## [1] "Run Completed at 2016-08-30 23:40:27"
#Load in data from Observed2m_Generate.Rmd
load("Observed.RData")
#load("ObservedModel.RData")
newModel<-T
#source functions
source("Bayesian/BayesFunctions.R")

View Raw Data

ggplot(indat,aes(x=Traitmatch,y=Yobs,col=as.factor(R))) + facet_wrap(~Hummingbird,ncol=4,scales="free") + geom_point() + geom_smooth(method = "glm",method.args=list(family="poisson")) + scale_color_manual("Resource Availability",values=c("black",'blue','red'))

ggplot(indat[indat$Yobs==1,],aes(x=Traitmatch,fill=as.factor(R))) + facet_wrap(~Hummingbird,ncol=4,scales="free") + geom_density(alpha=.7) + scale_fill_manual("Resource Availability",values=c("black",'blue','red'))

m<-glmer(data=indat[,],Yobs~Traitmatch*R+(1|Hummingbird),family="poisson")
summary(m)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Yobs ~ Traitmatch * R + (1 | Hummingbird)
##    Data: indat[, ]
## 
##      AIC      BIC   logLik deviance df.resid 
##   7146.0   7191.3  -3566.0   7132.0     4756 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4850 -0.5293 -0.3050 -0.1925 28.1772 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Hummingbird (Intercept) 1.656    1.287   
## Number of obs: 4763, groups:  Hummingbird, 19
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -1.079336   0.306095  -3.526 0.000422 ***
## Traitmatch         -0.033845   0.004302  -7.867 3.62e-15 ***
## RMedium            -0.261394   0.094822  -2.757 0.005839 ** 
## RHigh              -0.303452   0.087808  -3.456 0.000549 ***
## Traitmatch:RMedium -0.002927   0.006711  -0.436 0.662739    
## Traitmatch:RHigh    0.001511   0.006168   0.245 0.806440    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmtc RMedim RHigh  Trt:RM
## Traitmatch  -0.154                            
## RMedium     -0.147  0.485                     
## RHigh       -0.161  0.533  0.509              
## Trtmtch:RMd  0.101 -0.643 -0.734 -0.347       
## Trtmtch:RHg  0.114 -0.711 -0.350 -0.735  0.462
indat$pred<-predict(m,type="response")
ggplot(indat,aes(x=Traitmatch)) + geom_point(aes(y=Yobs)) + geom_line(aes(y=pred,col=as.factor(R)))

Hierarchical Model

For hummingbird species i feeding on plant species j observed at time k and sampling event observed by transect (YTransect) or camera (YCamera)

Observation Model:

\[ Yobs_{i,j,k,d} \sim Binomial(N_{i,j,k},\phi) \]

Process Model:

\[ N_{i,j,k} \sim Poisson(\lambda_{i,j,k}) \] \[ log(\lambda_{i,j,k}) = \alpha_i + \beta_{1,i} * Traitmatch_{i,j} + \beta_{2,i} *Resources_k + \beta_{3,i} * Traitmatch_{i,j} * Resources_k \]

Priors

\[ \phi \sim Normal(0,0.386) \] \[\alpha_i \sim Normal(\alpha_\mu,\alpha_\tau)\] \[\beta_{1,i} \sim Normal(\mu_{\beta_1,\tau_{beta_1}})\] \[\beta_{2,i} \sim Normal(\mu_{\beta_2,\tau_{beta_2}})\] \[\beta_{3,i} \sim Normal(\mu_{\beta_3,\tau_{beta_3}})\]

Group Level Logit Transformed Means \[ \mu_\alpha \sim Normal(0,1.67)\] \[\mu_{\beta_1} \sim Normal(0,1.67)\] \[\mu_{\beta_2} \sim Normal(0,1.67)\] \[\mu_{\beta_3} \sim Normal(0,1.67)\]

Group Level Variance \[\tau_{\alpha} \sim Half cauchy(0,1,1)\] \[\tau_{\beta_1} \sim Half cauchy(0,1,1)\] \[\tau_{\beta_2} \sim Half cauchy(0,1,1)\] \[\tau_{\beta_3} \sim Half cauchy(0,1,1)\]

#Source model
source("Bayesian/NmixturePoissonRagged2m.R")

#print model
writeLines(readLines("Bayesian/NmixturePoissonRagged2m.R"))
## 
## sink("Bayesian/NmixturePoissonRagged2m.jags")
## 
## cat("
##     model {
##     #Compute true state for each pair of birds and plants
##     for (i in 1:Birds){
##     for (j in 1:Plants){
##     for (k in 1:Times){
##     
##     #Process Model
##     log(rho[i,j,k])<-alpha[i] + beta1[i] * Traitmatch[i,j] + beta2[i] * resources[i,j,k] + beta3[i] * resources[i,j,k] * Traitmatch[i,j] 
##     
##     #True State
##     S[i,j,k] ~ dpois(rho[i,j,k])
##     }
##     }
##     }
##     
##     #Observation Model
##     for (x in 1:Nobs){
##     
##     #Observation Process for cameras
##     Yobs_camera[x] ~ dbin(detect[Bird[x]],S[Bird[x],Plant[x],Time[x]])    
## 
##     #     #Assess Model Fit - Posterior predictive check
##     # 
##     #     #Fit discrepancy statistics
##          #Camera
##          eval_cam[x]<-detect[Bird[x]]*S[Bird[x],Plant[x],Time[x]]
##          E[x]<-pow((Yobs_camera[x]-eval_cam[x]),2)/(eval_cam[x]+0.5)
##     #     
##          ynew_cam[x]~dbin(detect[Bird[x]], S[Bird[x],Plant[x],Time[x]])
##          E.new[x]<-pow((ynew_cam[x]-eval_cam[x]),2)/(eval_cam[x]+0.5)
##     }
##     
##     #Priors
##     #Observation model
##     #Detect priors, logit transformed - Following lunn 2012 p85
##     
##     for(x in 1:Birds){
##     #For Cameras
##     logit(detect[x])<-dcam[x]
##     dcam[x]~dnorm(omega_mu,omega_tau)
##     }
##     
##     
##     #Process Model
##     #Species level priors
##     for (i in 1:Birds){
##     
##     #Intercept
##     #logit prior, then transform for plotting
##     alpha[i] ~ dnorm(alpha_mu,alpha_tau)
##     
##     #Traits slope 
##     beta1[i] ~ dnorm(beta1_mu,beta1_tau)    
##     
##     #Plant slope
##     beta2[i] ~ dnorm(beta2_mu,beta2_tau)    
##     
##     #Interaction slope
##     beta3[i] ~ dnorm(beta3_mu,beta3_tau)
## 
##     }
##     
##     #OBSERVATION PRIOR
##     omega_mu ~ dnorm(0,0.386)
##     omega_tau ~ dunif(0,10)
##     
##     #Group process priors
##     
##     #Intercept 
##     alpha_mu ~ dnorm(0,0.386)
##     alpha_tau ~ dt(0,1,1)I(0,)
##     alpha_sigma<-pow(1/alpha_tau,0.5) 
##     
##     #Trait
##     beta1_mu~dnorm(0,0.386)
##     beta1_tau ~ dt(0,1,1)I(0,)
##     beta1_sigma<-pow(1/beta1_tau,0.5)
##     
##     #Resources
##     beta2_mu~dnorm(0,0.386)
##     beta2_tau ~ dt(0,1,1)I(0,)
##     beta2_sigma<-pow(1/beta2_tau,0.5)
##     
##     #Interaction
##     beta3_mu~dnorm(0,0.386)
##     beta3_tau ~ dt(0,1,1)I(0,)
##     beta3_sigma<-pow(1/beta3_tau,0.5)
## 
##     #derived posterior check
##     fit<-sum(E[]) #Discrepancy for the observed data
##     fitnew<-sum(E.new[])
##     
##     }
##     ",fill=TRUE)
## 
## sink()
#Inits
InitStage <- function(){
  #A blank Y matrix - all present
  initY<-array(dim=c(Birds,Plants,Times),22)
  initB<-rep(0.5,Birds)
list(S=initY,dcam=initB)}

#Parameters to track
ParsStage <- c("alpha","beta1","beta2","beta3","alpha_mu","alpha_sigma","beta1_sigma","beta1_mu","beta2_mu","beta2_sigma","beta3_mu","beta3_sigma","detect","E","fit","fitnew")


#Jags Data
Dat<-list(
  Yobs_camera = indat$Camera,
  Birds=max(indat$jBird),
  Bird=indat$jBird,
  Plant=indat$jPlant,
  Time=indat$jID,
  Plants=max(indat$jPlant),
  Times=max(indat$jID),
  resources=resourceMatrix,
  Nobs=nrow(indat),
  Traitmatch=jTraitmatch)

  #MCMC options
  if(newModel){
    system.time(
      m2<-jags.parallel(data=Dat,parameters.to.save =ParsStage,inits=InitStage,model.file="Bayesian/NmixturePoissonRagged2m.jags",n.thin=10,n.iter=300000,n.burnin=295000,n.chains=2,DIC=F)
      )
  }
##      user    system   elapsed 
##     9.944     0.359 14286.424
#recompile if needed
runs<-200000

recompile(m2)

if(!newModel){
  system.time(m2<-update(m2,n.iter=runs,n.burnin=runs-5000,n.thin=10))  
}
#extract par to data.frame
pars_detect<-extract_par(m2,data=indat,Bird="jBird",Plant="jPlant")

Assess Convergence

###Chains
ggplot(pars_detect[pars_detect$par %in% c("alpha","beta1","beta2","beta3"),],aes(x=Draw,y=estimate,col=as.factor(Chain))) + geom_line() + facet_grid(par~species,scale="free") + theme_bw() + labs(col="Chain") + ggtitle("Species Level")

ggplot(pars_detect[pars_detect$par %in% c("detect"),],aes(x=Draw,y=estimate,col=as.factor(Chain))) + geom_line() + facet_wrap(~species,scale="free",ncol=3) + theme_bw() + labs(col="Chain") + ggtitle("Species Level")

ggplot(pars_detect[pars_detect$par %in% c("beta1_mu","beta1_sigma","beta2_mu","beta2_sigma","beta3_mu","beta3_sigma","alpha_mu","alpha_sigma"),],aes(x=Draw,y=estimate,col=as.factor(Chain))) + geom_line() + theme_bw() + labs(col="Chain") + ggtitle("Group Level Parameters") + facet_wrap(~par,scales="free")

Posteriors

###Posterior Distributions
ggplot(pars_detect[pars_detect$par %in% c("alpha","beta1","beta2","beta3"),],aes(x=estimate)) + geom_histogram(position='identity') +  facet_grid(species~par,scales="free") + theme_bw() + ggtitle("Species Parameters")

#Detection figure
pars_detect<-merge(pars_detect,jagsIndexBird,by.x="species",by.y="jBird",all=T)
ggplot(pars_detect[pars_detect$par %in% "detect",],aes(x=par,y=estimate)) + geom_violin(fill='black') + theme_bw() + ggtitle("Detection Probability") + facet_wrap(~Hummingbird)

Detection Table

detecttable<-group_by(pars_detect,Hummingbird,par) %>% filter(par %in% c('detect')) %>% summarize(mean=mean(estimate),lower=quantile(estimate,0.05),upper=quantile(estimate,0.95))
detecttable
## Source: local data frame [19 x 5]
## Groups: Hummingbird [?]
## 
##                  Hummingbird    par      mean      lower     upper
##                       (fctr) (fctr)     (dbl)      (dbl)     (dbl)
## 1             Andean Emerald detect 0.2460116 0.13384758 0.3994516
## 2         Booted Racket-tail detect 0.2204931 0.13802587 0.3246381
## 3                 Brown Inca detect 0.1666781 0.10561342 0.2510915
## 4        Buff-tailed Coronet detect 0.1877761 0.10501788 0.2843066
## 5              Collared Inca detect 0.1970451 0.10464681 0.3110083
## 6          Crowned Woodnymph detect 0.2204191 0.11881498 0.3418982
## 7    Fawn-breasted Brilliant detect 0.1928675 0.08611181 0.3428303
## 8          Gorgeted Sunangel detect 0.4460895 0.27255985 0.6511495
## 9    Green-crowned Brilliant detect 0.2229545 0.08600007 0.4159489
## 10   Green-fronted Lancebill detect 0.2741602 0.13855144 0.4100603
## 11             Hoary Puffleg detect 0.2360042 0.11421216 0.3988843
## 12    Purple-bibbed Whitetip detect 0.2396576 0.10982839 0.4048972
## 13 Rufous-tailed Hummingbird detect 0.2073415 0.08644013 0.3740578
## 14      Speckled Hummingbird detect 0.1882605 0.10919659 0.2823925
## 15    Stripe-throated Hermit detect 0.2368006 0.13394028 0.3496642
## 16      Tawny-bellied Hermit detect 0.3198216 0.19842256 0.4544295
## 17       Violet-tailed Sylph detect 0.2240017 0.11694238 0.3533654
## 18  Wedge-billed Hummingbird detect 0.1568365 0.07167055 0.2771828
## 19    White-whiskered Hermit detect 0.2468784 0.15660755 0.3671323

Discrepancy

The goodness of fit is a measured as chi-squared. The expected value for each day is the detection rate * the estimate intensity of interactions. The expected value is compared to the observed value of the actual data. In addition, a replicate dataset is generated from the posterior predicted intensity. Better fitting models will have lower discrepancy values and be Better fitting models are smaller values and closer to the 1:1 line. A perfect model would be 0 discrepancy. This is unrealsitic given the stochasticity in the sampling processes. Rather, its better to focus on relative discrepancy. In addition, a model with 0 discrepancy would likely be seriously overfit and have little to no predictive power.

fitstat<-pars_detect[pars_detect$par %in% c("fit","fitnew"),]
fitstat<-dcast(fitstat,Draw+Chain~par,value.var="estimate")

ymin<-round(min(fitstat$fit))
ymax<-round(max(fitstat$fit))
ab<-data.frame(x=0:ymax,y=0:ymax)
disc_obs<-ggplot(fitstat,aes(x=fit,y=fitnew)) + geom_point() + theme_bw() + labs(x="Discrepancy of observed data",y="Discrepancy of replicated data",col="Model")  + ggtitle("Posterior predictive check") + geom_line(data=ab,aes(x=x,y=y)) + coord_fixed() + ylim(ymin=0,ymax=max(max(c(fitstat$fit,fitstat$fitnew)))) + xlim(xmin=0,xmax=max(max(c(fitstat$fit,fitstat$fitnew))))
disc_obs

#Bayesian p-value
sum(fitstat$fitnew>fitstat$fit)/nrow(fitstat)
## [1] 0
ggsave("Figures/ObservedDiscrepancy.jpeg",width = 5,height=10)

Predicted Relationship

#Expand out pars
castdf<-dcast(pars_detect[pars_detect$par %in% c("beta1_mu","beta2_mu","beta3_mu","alpha_mu"),], Chain + Draw~par,value.var="estimate")

Posterior prediction

#Trajectories from posterior
predy<-trajF(alpha=castdf$alpha_mu,beta1=castdf$beta1_mu,trait=indat$Traitmatch,resources=indat$scaledR,beta2=castdf$beta2_mu,beta3=castdf$beta3_mu)

ggplot(data=predy,aes(x=trait)) + geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.4,fill="red")  +  theme_bw() + ylab("Interactions") + xlab("Difference between Bill and Corolla Length") + geom_point(data=indat,aes(x=Traitmatch,y=Camera)) + geom_line(aes(y=mean)) + geom_point(data=indat,aes(x=Traitmatch,y=Transect)) 

At High and Low Resource Availability

#Trajectories from posterior
predH<-trajF(alpha=castdf$alpha_mu,beta1=castdf$beta1_mu,trait=indat[indat$R=="High","Traitmatch"],resources=indat[indat$R=="High","scaledR"],beta2=castdf$beta2_mu,beta3=castdf$beta3_mu)

predM<-trajF(alpha=castdf$alpha_mu,beta1=castdf$beta1_mu,trait=indat[indat$R=="Medium","Traitmatch"],resources=indat[indat$R=="Medium","scaledR"],beta2=castdf$beta2_mu,beta3=castdf$beta3_mu)

predL<-trajF(alpha=castdf$alpha_mu,beta1=castdf$beta1_mu,trait=indat[indat$R=="Low","Traitmatch"],resources=indat[indat$R=="Low","scaledR"],beta2=castdf$beta2_mu,beta3=castdf$beta3_mu)

predhl<-melt(list(High=predH,Medium=predM,Low=predL),id.vars=colnames(predH))

colnames(predhl)[6]<-"BFlowerL"

predhl$BFlowerL<-factor(as.character(predhl$BFlowerL))
levels(predhl$BFlowerL)<-c("Low","Medium","High")

ggplot(data=predhl,aes(x=trait)) + geom_ribbon(aes(ymin=lower,ymax=upper,fill=BFlowerL),alpha=0.5)  + geom_line(aes(y=mean,col=BFlowerL),size=.8) + theme_bw() + ylab("Interactions") + xlab("Difference between Bill and Corolla Length") + geom_point(data=mindat,aes(x=Traitmatch,y=value))+ labs(fill="Resource Availability",col="Resource Availability") 

ggsave("Figures/AllRegression.jpeg",height=5,width=7)

Species Predictions

castdf<-dcast(pars_detect[pars_detect$par %in% c("beta1","beta2","beta3","alpha"),], species +Chain +Draw ~par ,value.var="estimate")

#Turn to 
castdf$species<-factor(castdf$species,levels=1:max(as.numeric(castdf$species)))

species.split<-split(castdf,list(castdf$species),drop = T)

species.traj<-list()

for(d in 1:length(species.split)){
  x<-species.split[[d]]
  index<-unique(x$species)
  
  #get data for those species
  billd<-indat[indat$jBird %in% index,]

  #scale resources
  species.traj[[d]]<-trajF(alpha=x$alpha,beta1=x$beta1,beta2=x$beta2,beta3=x$beta3,resources=billd$scaledR,trait=billd$Traitmatch)
  }

names(species.traj)<-names(species.split)

species.traj<-melt(species.traj,id.var=colnames(species.traj[[1]]))

#split out names and model
species.traj[,c("Index")]<-colsplit(species.traj$L1,"\\.",c("Index"))

spe<-merge(species.traj,jagsIndexBird,by.x="Index",by.y="jBird")

#plot and compare to original data
ggplot(data=spe,aes(x=trait)) + geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.2,fill='red') + geom_line(aes(y=mean),size=.5) + theme_bw() + ylab("Daily Interaction Rate")+ xlab("Difference between Bill and Corolla Length") + facet_wrap(~Hummingbird,scales="free",ncol=4) + geom_point(data=mindat,aes(x=Traitmatch,y=value),size=2.5) 

Species Predictions: High and Low Availability

castdf<-dcast(pars_detect[pars_detect$par %in% c("beta1","beta2","beta3","alpha"),], species +Chain +Draw ~par ,value.var="estimate")

#Turn to 
castdf$species<-factor(castdf$species,levels=1:max(as.numeric(castdf$species)))

species.split<-split(castdf,list(castdf$species),drop = T)

species.traj<-list()

for(d in 1:length(species.split)){
  x<-species.split[[d]]
  index<-unique(x$species)
  
  #get data for those species
  billd<-indat[indat$jBird %in% index,]

  sh<-trajF(alpha=x$alpha,beta1=x$beta1,beta2=x$beta2,beta3=x$beta3,resources=billd[billd$R=="High","scaledR"],trait=billd[billd$R=="High","Traitmatch"])
  
  sm<-trajF(alpha=x$alpha,beta1=x$beta1,beta2=x$beta2,beta3=x$beta3,resources=billd[billd$R=="Medium","scaledR"],trait=billd[billd$R=="Medium","Traitmatch"])

  sl<-trajF(alpha=x$alpha,beta1=x$beta1,beta2=x$beta2,beta3=x$beta3,resources=billd[billd$R=="Low","scaledR"],trait=billd[billd$R=="Low","Traitmatch"])
  
  sm<-melt(list(High=sh,Medium=sm,Low=sl),id.vars=colnames(sl))
  species.traj[[d]]<-sm
  }

names(species.traj)<-names(species.split)

species.traj<-melt(species.traj,id.var=colnames(species.traj[[1]]))

#split out names and model
species.traj[,c("Index")]<-colsplit(species.traj$L1,"\\.",c("Index"))

spe<-merge(species.traj,jagsIndexBird,by.x="Index",by.y="jBird")

#plot and compare to original data
indat$rplot<-as.factor(indat$scaledR)
levels(indat$rplot)<-c("Low","Medium","High")
spe$resources<-as.factor(spe$resources)
levels(spe$resources)<-c("Low","Medium","High")
ggplot(data=spe,aes(x=trait)) + geom_ribbon(aes(ymin=lower,ymax=upper,fill=resources),alpha=0.5) + geom_line(aes(y=mean,group=resources),size=.5) + theme_bw() + ylab("Daily Interaction Rate")+ xlab("Difference between Bill and Corolla Length") + facet_wrap(~Hummingbird,scales="free",ncol=4) + geom_point(data=indat,aes(x=Traitmatch,y=Yobs),size=1.2) + scale_fill_manual("Resource Availability",values=c("Grey70","Grey40","Black"))

ggsave("Figures/SpeciesRegression.jpeg",height=7,width=11)

Species Level Interaction

castdf<-dcast(pars_detect[pars_detect$par %in% c("beta1","beta2","beta3","alpha"),], species +Chain + Draw~par,value.var="estimate")

#Turn to 
castdf$species<-factor(castdf$species,levels=1:max(as.numeric(castdf$species)))

species.split<-split(castdf,list(castdf$species),drop = T)

species.traj<-list()

for(d in 1:length(species.split)){
  dat<-species.split[[d]]
  index<-unique(dat$species)
  
  #get data for those species
  billd<-mindat[mindat$jBird %in% index,]

  #Calculate interaction effect
  species.traj[[d]]<-intF(alpha=dat$alpha,beta1=dat$beta1,x=billd[billd$value > 0 & !is.na(billd$value),'Traitmatch'],resources=billd[billd$value > 0 & !is.na(billd$value),'scaledR'],beta2=dat$beta2,beta3=dat$beta3)
  }

names(species.traj)<-names(species.split)
species.traj<-melt(species.traj,id.var=colnames(species.traj[[1]]))

#split out names and model
species.traj[,c("Index")]<-colsplit(species.traj$L1,"\\.",c("Index"))

spe<-merge(species.traj,jagsIndexBird,by.x="Index",by.y="jBird")

#match colnames

#plot and compare to original data
ggplot(data=spe,aes(x=x)) + geom_ribbon(aes(ymin=lower,ymax=upper,fill=Hummingbird),alpha=0.3) + geom_line(aes(y=mean,col=Hummingbird),size=1) + theme_bw() + xlab("Difference between Bill and Corolla Length")  + ylab("Effect of Resources on Trait Difference") + facet_wrap(~Hummingbird,scales="free",ncol=3)
ggsave("Figures/SpeciesInteraction.jpeg",height=6,width=7)

These plots can be tricky to interpret if one forgets that trait matching as a covariate is a distance. Therefore remember that a positive slope in the plot above indiciates, “As resources increase species use flowers less similiar to their bill lengths”.

Which species did we predict well?

smat<-pars_detect %>% filter(par=="E") %>% group_by(species) %>% summarize(fit=sum(estimate)) %>% arrange(desc(fit))
smat<-merge(smat,jagsIndexBird,by.x="species",by.y="jBird")

Discrepancy matrix

dmat<-pars_detect %>% filter(par=="E") %>% group_by(species,plant) %>% summarize(fit=sum(estimate)) %>% arrange(desc(fit))
dmat<-merge(dmat,jagsIndexBird,by.x="species",by.y="jBird")
dmat<-merge(dmat,jagsIndexPlants,by.x="plant",by.y="jPlant")
head(dmat,10)
##    plant species        fit             Hummingbird
## 1      1       4   28.50690     Buff-tailed Coronet
## 2      1       9   35.66932 Green-crowned Brilliant
## 3      1      17 2157.04299     Violet-tailed Sylph
## 4      1       7  921.59691 Fawn-breasted Brilliant
## 5      1      12   30.71657  Purple-bibbed Whitetip
## 6      1       6   60.13242       Crowned Woodnymph
## 7      1       3 8471.06069              Brown Inca
## 8      1       5  663.77643           Collared Inca
## 9      1      10  235.85366 Green-fronted Lancebill
## 10     1       2   69.23956      Booted Racket-tail
##                Iplant_Double
## 1  Alloplectus tetragonoides
## 2  Alloplectus tetragonoides
## 3  Alloplectus tetragonoides
## 4  Alloplectus tetragonoides
## 5  Alloplectus tetragonoides
## 6  Alloplectus tetragonoides
## 7  Alloplectus tetragonoides
## 8  Alloplectus tetragonoides
## 9  Alloplectus tetragonoides
## 10 Alloplectus tetragonoides
#get order
hord<-dmat %>% group_by(Hummingbird) %>% summarize(n=sum(fit)) %>% arrange(desc(n)) %>% .$Hummingbird
pord<-dmat %>% group_by(Iplant_Double) %>% summarize(n=sum(fit)) %>% arrange(n) %>%  .$Iplant_Double
dmat$Hummingbird<-factor(dmat$Hummingbird,levels=hord)
dmat$Iplant_Double<-factor(dmat$Iplant_Double,levels=pord)
ggplot(dmat,aes(x=Hummingbird,Iplant_Double,fill=fit)) + geom_tile() + theme_bw() + scale_fill_continuous("Discrepancy",low="blue",high="red") + theme(axis.text.x=element_text(angle = -90))

Resource Abundance Functions

Let’s take a closer look at distribution of interaction effect posteriors values for each species.

post<-pars_detect %>% filter(par %in% "beta2") %>% group_by(Hummingbird) %>% summarize(mean=mean(estimate),median=median(estimate),lower=quantile(probs=0.025,estimate),upper=quantile(probs=0.975,estimate)) %>% melt(id.vars='Hummingbird')
ggplot(pars_detect[pars_detect$par %in% "beta2",],aes(x=estimate)) + geom_histogram() + facet_wrap(~Hummingbird,scales='free',ncol=4) + geom_vline(data=post,aes(xintercept=value,col=variable))

Interaction density functions

Let’s take a closer look at distribution of interaction effect posteriors values for each species.

post<-pars_detect %>% filter(par %in% "beta3") %>% group_by(Hummingbird) %>% summarize(mean=mean(estimate),median=median(estimate),lower=quantile(probs=0.025,estimate),upper=quantile(probs=0.975,estimate)) %>% melt(id.vars='Hummingbird')
ggplot(pars_detect[pars_detect$par %in% "beta3",],aes(x=estimate)) + geom_histogram() + facet_wrap(~Hummingbird,scales='free',ncol=4) + geom_vline(data=post,aes(xintercept=value,col=variable))

Trait-matching and Bill Length

Do species with long bill lengths have positive traitmatching effects?

#species names
b<-pars_detect[pars_detect$par %in% "beta1",]

#traits
b<-merge(b,hum.morph,by.x="Hummingbird",by.y="English")

post<-b %>% filter(par %in% "beta1") %>% group_by(Hummingbird) %>% summarize(mean=mean(estimate),median=median(estimate),lower=quantile(probs=0.025,estimate),upper=quantile(probs=0.975,estimate),quantile_l=quantile(estimate)[[1]],quantile_u=quantile(estimate)[[2]]) %>% melt(id.vars='Hummingbird')

#get order of mean posterior
ord<-post %>% filter(variable=="mean") %>% arrange(value) %>% .$Hummingbird

b$Hummingbird<-factor(b$Hummingbird,levels=ord)
ggplot(b,aes(y=estimate,x=as.factor(Total_Culmen),group=Hummingbird)) + geom_violin(fill='grey50') + coord_flip()  + ggtitle("Trait-matching and Bill Length") + theme_bw()

Interaction and Bill Length

Do species with long bill lengths have positive interaction effects?

#species names
b<-pars_detect[pars_detect$par %in% "beta3",]

#traits
b<-merge(b,hum.morph,by.x="Hummingbird",by.y="English")

post<-b %>% filter(par %in% "beta3") %>% group_by(Hummingbird) %>% summarize(mean=mean(estimate),median=median(estimate),lower=quantile(probs=0.025,estimate),upper=quantile(probs=0.975,estimate),quantile_l=quantile(estimate)[[1]],quantile_u=quantile(estimate)[[2]]) %>% melt(id.vars='Hummingbird')

#get order of mean posterior
ord<-post %>% filter(variable=="mean") %>% arrange(value) %>% .$Hummingbird

b$Hummingbird<-factor(b$Hummingbird,levels=ord)
ggplot(b,aes(y=estimate,x=Hummingbird,fill=Total_Culmen)) + geom_violin() + coord_flip() + scale_fill_continuous(low='blue',high='red') + ggtitle("Interaction Effect and Bill Length") + theme_bw()

Estimated niche breadth

castdf<-dcast(pars_detect[pars_detect$par %in% c("beta1","beta2","beta3","alpha"),], species +Chain + Draw~par,value.var="estimate")

#Turn to 
castdf$species<-factor(castdf$species,levels=1:max(as.numeric(castdf$species)))

species.split<-split(castdf,list(castdf$species),drop = T)

species.traj<-lapply(species.split,function(dat){
  index<-unique(dat$species)
  
  #get data for those species
  billd<-indat[indat$jBird %in% index,]
  
  d<-data.frame(alpha=dat$alpha,beta1=dat$beta1,beta2=dat$beta2,beta3=dat$beta3)
  
  #fit regression for each input estimate
  sampletraj<-list()
  
  for (y in 1:nrow(d)){
    v=exp(d$alpha[y] + d$beta1[y] * billd$Traitmatch + d$beta2[y] * billd$scaledR + d$beta3[y] * billd$Traitmatch*billd$scaledR)
    
    sampletraj[[y]]<-data.frame(x=as.numeric(billd$Traitmatch),y=as.numeric(v),r=as.numeric(billd$scaledR),jBird=billd$jBird,jPlant=billd$jPlant,jTime=billd$jTime)
  }
  
  sample_all<-rbind_all(sampletraj)
})
  
species.traj<-rbind_all(species.traj)

Mean Estimates for Corolla Sizes

species.mean<-species.traj %>% group_by(jBird,jPlant,r) %>% summarize(Traitmatch=unique(x),phi=mean(y))

tomerge<-indat %>% select(jBird,jPlant,Hummingbird,Iplant_Double) %>% distinct()

species.mean<-merge(species.mean,tomerge)

#get corolla sizes
species.mean<-merge(species.mean,fl.morph,by.x="Iplant_Double", by.y="Group.1")

#bill order
ord<-hum.morph %>% arrange(Total_Culmen) %>% .$English
species.mean$Hummingbird<-factor(species.mean$Hummingbird,levels=ord)

#add level to hum.morph to match naming convention
species.mean<-merge(species.mean,hum.morph[,c("English","Total_Culmen")],by.x="Hummingbird",by.y="English")

p<-ggplot(species.mean,aes(x=Corolla,y=phi,col=as.factor(r))) + geom_line(size=.9) + geom_vline(aes(xintercept=Total_Culmen),linetype='dashed') + facet_wrap(~Hummingbird,ncol=3,scales="free_y")  + theme_bw() + ylab("Probability of Interaction") + scale_color_manual("Resource Availability",labels=c("Low","Medium","High"),values=c("Black","Blue","Red")) + xlab("Flower Corolla Length (mm)") 
p

ggsave("Figures/ResponseCurves.jpeg",height=6.5,width=9)

Niche Breadth

species.mean<-species.traj %>% group_by(jBird,jPlant,r) %>% summarize(Traitmatch=unique(x),phi=mean(y),phi_low=quantile(y,0.05),phi_high=quantile(y,0.95))

#merge names
species.mean<-merge(species.mean,jagsIndexBird)
species.mean<-merge(species.mean,jagsIndexPlants)

#get corolla sizes
species.mean<-merge(species.mean,fl.morph,by.x="Iplant_Double", by.y="Group.1")

#bill order
ord<-hum.morph %>% arrange(Total_Culmen) %>% .$English
species.mean$Hummingbird<-factor(species.mean$Hummingbird,levels=ord)

#add level to hum.morph to match naming convention
species.mean<-merge(species.mean,hum.morph[,c("English","Total_Culmen")],by.x="Hummingbird",by.y="English")

#label factor
species.mean$r<-as.factor(species.mean$r)
levels(species.mean$r)<-c("Low","Medium","High")
ggplot(species.mean) + geom_ribbon(alpha=0.7,aes(x=Corolla,ymin=phi_low,ymax=phi_high,fill=r)) + theme_bw() + facet_wrap(~Hummingbird,scales="free",ncol=4)+ ggtitle("Niche Breadth") + geom_vline(aes(xintercept=Total_Culmen),linetype='dashed') + geom_line(aes(x=Corolla,y=phi,fill=r)) + ylab("Probability of Interaction") + xlab("Corolla Length (mm)") + scale_fill_manual("Resource Availability",values=c("Grey70","Grey20","Black"))

ggsave("Figures/NicheBreadth.jpeg",height=7,width=11)

Range plots

nsplit<-split(species.mean,species.mean$r)
makeR<-function(x){
  
  #input matrix
  aggm<-matrix(nrow=nrow(jagsIndexBird),ncol=nrow(jagsIndexPlants),data=0)
  for (j in 1:nrow(x)){
    aggm[x[j,"jBird"],x[j,"jPlant"]]<-rpois(1,lambda=x[j,"phi"])
  }
  
  aggm<-melt(aggm)
  colnames(aggm)<-c("jBird","jPlant","P")
  tomerge<-species.mean %>% select(jBird,jPlant,Corolla,Hummingbird,Traitmatch) %>% distinct()
  aggm<-merge(aggm,tomerge)
  return(aggm)
}

nstat<-lapply(nsplit,function(x){
  netstat<-melt(lapply(1:100,function(k) makeR(x)),id.vars=c("jBird","jPlant","Hummingbird","Traitmatch","Corolla","P"))
  colnames(netstat)[7]<-"Iteration"
  return(netstat)
})

names(nstat)<-c("Low","Medium","High")
nstat<-melt(nstat,colnames(nstat[[1]]))

Predicted density

#order the levels
nstat$L1<-factor(nstat$L1,levels=c("High","Medium","Low"))

#add level to hum.morph to match naming convention
nstat<-merge(nstat,hum.morph[,c("English","Total_Culmen")],by.x="Hummingbird",by.y="English")

ggplot(nstat[nstat$P>0,],aes(x=Corolla,fill=L1)) + geom_density(alpha=0.8) + facet_wrap(~Hummingbird,scales='free',ncol=4)  + theme_bw()+ scale_fill_manual("Resource Availability",values=c("Black","Grey40","Grey85"))  + geom_vline(aes(xintercept=Total_Culmen),linetype='dashed')

ggsave("Figures/PredictedDensity.jpeg",dpi=400,height=8,width=10)
rangestat<-nstat %>% filter(P==1) %>% group_by(Hummingbird,L1) %>% summarize(mean=mean(Corolla),lower=quantile(Corolla,0.05),upper=quantile(Corolla,0.95))
  
ggplot(rangestat,aes(x=Hummingbird,ymin=lower,ymax=upper,ymean=mean,col=L1)) + geom_linerange(alpha=0.5,position=position_dodge(width=.75),size=2) + geom_point(aes(y=mean),position=position_dodge(width=.75)) + labs(col="Resource Availability",y="Corolla Length (mm)") + theme_bw() + coord_flip() + theme(legend.position="bottom")

Predicted Mean

meanstat<-nstat %>% filter(P==1) %>% group_by(Hummingbird,L1,Iteration) %>% summarize(a=mean(Corolla))%>% summarize(mean=mean(a),lower=quantile(a,0.05),upper=quantile(a,0.95))
  
ggplot(meanstat,aes(x=Hummingbird,ymin=lower,ymax=upper,ymean=mean,col=L1)) + geom_linerange(alpha=0.5,position=position_dodge(width=.6),size=2) + geom_point(aes(y=mean),position=position_dodge(width=.6)) + labs(col="Resource Availability",y="Corolla Length (mm)") + theme_bw() + coord_flip() + theme(legend.position="bottom")

Predicted Min

meanstat<-nstat %>% filter(P==1) %>% group_by(Hummingbird,L1,Iteration) %>% summarize(a=min(Corolla))%>% summarize(mean=mean(a),lower=quantile(a,0.05),upper=quantile(a,0.95))
  
ggplot(meanstat,aes(x=Hummingbird,ymin=lower,ymax=upper,ymean=mean,col=L1)) + geom_linerange(alpha=0.5,position=position_dodge(width=.75),size=2) + geom_point(aes(y=mean),position=position_dodge(width=.75)) + labs(col="Resource Availability",y="Corolla Length (mm)") + theme_bw() + coord_flip() + theme(legend.position="bottom")

Predicted Max

meanstat<-nstat %>% filter(P==1) %>% group_by(Hummingbird,L1,Iteration) %>% summarize(a=min(Corolla))%>% summarize(mean=mean(a),lower=quantile(a,0.05),upper=quantile(a,0.95))
  
ggplot(meanstat,aes(x=Hummingbird,ymin=lower,ymax=upper,ymean=mean,col=L1)) + geom_linerange(alpha=0.5,position=position_dodge(width=.6),size=2) + geom_point(aes(y=mean),position=position_dodge(width=.6)) + labs(col="Resource Availability",y="Corolla Length (mm)") + theme_bw() + coord_flip() + theme(legend.position="bottom")

Generate network

#Split by resource
nsplit<-split(species.mean,species.mean$r)

makeN<-function(x){
  
  #input matrix
  aggm<-matrix(nrow=nrow(jagsIndexBird),ncol=nrow(jagsIndexPlants),data=0)
  for (j in 1:nrow(x)){
    aggm[x[j,"jBird"],x[j,"jPlant"]]<-rpois(1,lambda=x[j,"phi"])
  }
  #calculate network statistic
  nstat<-networklevel(aggm,index=c("niche overlap"),level="lower")
}

nstat<-lapply(nsplit,function(x){
  netstat<-melt(t(sapply(1:500,function(k) makeN(x)))) 
  colnames(netstat)<-c("Iteration","Metric","value")
  return(netstat)
})

names(nstat)<-c("Low","Medium","High")
nstat<-melt(nstat,colnames(nstat[[1]]))

nstat$L1<-factor(nstat$L1,c("Low","Medium","High"))
levels(nstat$Metric)<-"Niche Overlap"

ggplot(nstat,aes(x=value,fill=L1)) + geom_density(alpha=0.7) + facet_wrap(~Metric,scales='free',nrow=3) + scale_fill_discrete("Resource Availability") + theme_bw() + scale_fill_manual("Resource Availability",values=c("Grey80","Grey50","Black"))

ggsave("Figures/NetworkStatistics.jpeg",height=4,width=6,dpi=600) 

Compared to raw visits

sindat<-split(indat,indat$scaledR)
ndat<-lapply(sindat,function(x){
  web<-acast(x[x$Yobs==1,],Hummingbird~Iplant_Double,value.var = "Yobs", fun= function(x){ length(unique(x))})
  nstat<-t(networklevel(web,index=c("niche overlap"),level="lower"))
})
names(ndat)<-c("Low","Medium","High")
ndat<-melt(ndat)

colnames(ndat)<-c("Var1","Metric","value","L1")


ggplot(nstat,aes(x=value,fill=L1)) + geom_density(alpha=0.6) + facet_wrap(~Metric,scales='free',nrow=2) + scale_fill_discrete("Resource Availability") + geom_vline(data=ndat,aes(xintercept=value,col=L1),linetype="dashed") + scale_color_discrete(guide='none') 

save.image("ObservedModel.RData")